library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.3
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.8
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.3
## Warning: package 'readr' was built under R version 4.1.3
## Warning: package 'dplyr' was built under R version 4.1.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(gtable)
library(gtsummary)
## Warning: package 'gtsummary' was built under R version 4.1.3
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.1.3
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(knitr)
library(car)
## Warning: package 'car' was built under R version 4.1.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.1.3
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
library(boot)
## Warning: package 'boot' was built under R version 4.1.3
##
## Attaching package: 'boot'
## The following object is masked from 'package:car':
##
## logit
data <- readRDS("Dataset.RData")
data$y = 0
data$y[data$survival == 'Long'] = 1
# immune selected genes MC1R, AKT1, NAMPT, TNFRSF19, PAK2, SOS1, B2M
l1 <- glm(data$y ~ data$MC1R + as.factor(data$stage), family = 'binomial')
l1$coefficients[2]
## data$MC1R
## -0.9632995
l2 <- glm(data$y ~ data$AKT1 + as.factor(data$stage), family = 'binomial')
l2$coefficients[2]
## data$AKT1
## -0.917543
l3 <- glm(data$y ~ data$NAMPT + as.factor(data$stage), family = 'binomial')
l3$coefficients[2]
## data$NAMPT
## 0.1939064
l4 <- glm(data$y ~ data$TNFRSF19 + as.factor(data$stage), family = 'binomial')
l4$coefficients[2]
## data$TNFRSF19
## 0.23159
l5 <- glm(data$y ~ data$PAK2 + as.factor(data$stage), family = 'binomial')
l5$coefficients[2]
## data$PAK2
## 0.5984475
l6 <- glm(data$y ~ data$SOS1 + as.factor(data$stage), family = 'binomial')
l6$coefficients[2]
## data$SOS1
## 0.2604371
l7 <- glm(data$y ~ data$B2M + as.factor(data$stage), family = 'binomial')
l7$coefficients[2]
## data$B2M
## 0.3445805
dd = data[c('MC1R', 'AKT1', 'NAMPT', 'TNFRSF19', 'PAK2', 'SOS1', 'B2M')]
library(plotly)
## Warning: package 'plotly' was built under R version 4.1.3
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(ggcorrplot)
## Warning: package 'ggcorrplot' was built under R version 4.1.3
library(heatmaply)
## Warning: package 'heatmaply' was built under R version 4.1.3
## Loading required package: viridis
## Loading required package: viridisLite
##
## ======================
## Welcome to heatmaply version 1.3.0
##
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
##
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## You may ask questions at stackoverflow, use the r and heatmaply tags:
## https://stackoverflow.com/questions/tagged/heatmaply
## ======================
heatmaply_cor(cor(dd))
cor(dd)
## MC1R AKT1 NAMPT TNFRSF19 PAK2 SOS1
## MC1R 1.0000000 0.6166103 -0.1633626 -0.13110913 -0.3211794 -0.35683222
## AKT1 0.6166103 1.0000000 -0.2861000 -0.24337733 -0.4490884 -0.37407132
## NAMPT -0.1633626 -0.2861000 1.0000000 0.23196342 0.1865976 0.11443004
## TNFRSF19 -0.1311091 -0.2433773 0.2319634 1.00000000 0.2679136 0.04096849
## PAK2 -0.3211794 -0.4490884 0.1865976 0.26791357 1.0000000 0.21673064
## SOS1 -0.3568322 -0.3740713 0.1144300 0.04096849 0.2167306 1.00000000
## B2M -0.2259642 -0.1928459 0.2008789 -0.08569906 0.1378465 0.16882559
## B2M
## MC1R -0.22596423
## AKT1 -0.19284591
## NAMPT 0.20087893
## TNFRSF19 -0.08569906
## PAK2 0.13784648
## SOS1 0.16882559
## B2M 1.00000000
dd1 = as.data.frame(dd)
dd1$y = data$y
dd1$age_c = data$age - mean(data$age)
dd1$sex = 0
dd1$sex[data$gender == 'female'] = 1
dd1$pathStage = as.factor(data$stage)
model = glm(y ~ MC1R + AKT1 + NAMPT + TNFRSF19 + PAK2 + SOS1 + B2M + pathStage, data = dd1, family = 'binomial')
vif(model)
## GVIF Df GVIF^(1/(2*Df))
## MC1R 1.199758 1 1.095335
## AKT1 1.435392 1 1.198078
## NAMPT 1.127966 1 1.062058
## TNFRSF19 1.183216 1 1.087757
## PAK2 1.197908 1 1.094490
## SOS1 1.192218 1 1.091888
## B2M 1.131922 1 1.063918
## pathStage 1.134084 3 1.021192
# selected top genes: MC1R, AKT1, PAK2, SOS1, B2M
# model 1: basline + selected top genes
model1 <- glm(y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + pathStage, data = dd1, family = 'binomial')
summary(model1)
##
## Call:
## glm(formula = y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c +
## pathStage, family = "binomial", data = dd1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5591 -0.9368 0.2876 0.9012 2.0907
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.83506 0.46974 3.907 9.36e-05 ***
## MC1R -0.51267 0.30252 -1.695 0.09014 .
## AKT1 -0.60824 0.24547 -2.478 0.01322 *
## PAK2 0.32723 0.18588 1.760 0.07834 .
## SOS1 -0.10986 0.17285 -0.636 0.52506
## B2M 0.22297 0.15674 1.423 0.15488
## sex -0.23477 0.31377 -0.748 0.45432
## age_c -0.02149 0.01077 -1.995 0.04600 *
## pathStage2 -1.64327 0.52583 -3.125 0.00178 **
## pathStage3 -2.33998 0.50959 -4.592 4.39e-06 ***
## pathStage4 -18.48491 882.74369 -0.021 0.98329
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 331.89 on 239 degrees of freedom
## Residual deviance: 251.47 on 229 degrees of freedom
## AIC: 273.47
##
## Number of Fisher Scoring iterations: 13
model1.diag <- glm.diag(model1)
data.frame("cooks" = model1.diag$cook, "leverage" = model1.diag$h) %>%
pivot_longer(everything(), names_to = "measure") %>%
mutate(cutoff = ifelse(measure=="cooks",1,26/240)) %>%
ggplot(aes(x = seq(0.5,240,by=0.5))) +
geom_line(aes(y = value, col = measure)) +
geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
theme_classic() +
facet_wrap(~measure, scales = "free_y", nrow = 2) +
scale_color_manual(values = c("orange","cornflowerblue")) +
labs(x = "observation")
p_res1 <- residuals(model1, type = "pearson")
d_res1 <- residuals(model1, type = "deviance")
data.frame("pearson" = p_res1, "deviance" = d_res1) %>%
pivot_longer(everything(), names_to = "residual") %>%
ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
geom_point() +
facet_wrap(~residual, nrow = 2) +
theme_classic() +
scale_color_manual(values = c("darkgreen","purple")) +
labs(x = "observation")
data.frame("pearson" = model1.diag$rp, "deviance" = model1.diag$rd) %>%
pivot_longer(everything(), names_to = "residual") %>%
ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
geom_point() +
facet_wrap(~residual, nrow = 2) +
theme_classic() +
scale_color_manual(values = c("darkgoldenrod2","orchid3")) +
labs(x = "observation", y = "standardized values")
#kable(data.frame("pearson" = p_res1, "deviance" = d_res1,"leverage" = model1.diag$h, "cooks" = model1.diag$cook,"p_standard" = model1.diag$rp, "d_standard" = model1.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
# model 2: baseline + selected top genes + immune resposne genes
dd1$CD86 = data$CD86
dd1$DDX58 = data$DDX58
dd1$KLRD1 = data$KLRD1
dd1$TLR4 = data$TLR4
dd1$DCK = data$DCK
dd1$C5 = data$C5
dd1$CCL28 = data$CCL28
dd1$RFX5 = data$RFX5
dd1$EIF2AK2 = data$EIF2AK2
model2 <- glm(y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4
+ C5 + CCL28 + RFX5 + EIF2AK2 + DCK + pathStage, data = dd1, family = 'binomial')
summary(model2)
##
## Call:
## glm(formula = y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c +
## CD86 + DDX58 + KLRD1 + TLR4 + C5 + CCL28 + RFX5 + EIF2AK2 +
## DCK + pathStage, family = "binomial", data = dd1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4960 -0.8346 0.2168 0.7899 2.1391
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.73513 0.48444 3.582 0.000341 ***
## MC1R -0.37369 0.32273 -1.158 0.246901
## AKT1 -0.41097 0.27016 -1.521 0.128214
## PAK2 0.22909 0.20232 1.132 0.257502
## SOS1 -0.26932 0.18330 -1.469 0.141755
## B2M -0.40808 0.26533 -1.538 0.124050
## sex -0.26646 0.33268 -0.801 0.423165
## age_c -0.03155 0.01228 -2.569 0.010209 *
## CD86 0.37462 0.32247 1.162 0.245341
## DDX58 -0.21035 0.22115 -0.951 0.341527
## KLRD1 0.14931 0.25890 0.577 0.564135
## TLR4 0.30899 0.24808 1.246 0.212944
## C5 0.19427 0.16974 1.145 0.252397
## CCL28 0.06700 0.19339 0.346 0.728989
## RFX5 0.30399 0.20599 1.476 0.140020
## EIF2AK2 0.32518 0.22064 1.474 0.140542
## DCK 0.46457 0.22625 2.053 0.040039 *
## pathStage2 -1.23248 0.56163 -2.194 0.028202 *
## pathStage3 -2.31474 0.53271 -4.345 1.39e-05 ***
## pathStage4 -17.61213 882.74403 -0.020 0.984082
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 331.89 on 239 degrees of freedom
## Residual deviance: 234.58 on 220 degrees of freedom
## AIC: 274.58
##
## Number of Fisher Scoring iterations: 13
model2.diag <- glm.diag(model2)
data.frame("cooks" = model2.diag$cook, "leverage" = model2.diag$h) %>%
pivot_longer(everything(), names_to = "measure") %>%
mutate(cutoff = ifelse(measure=="cooks",1,26/240)) %>%
ggplot(aes(x = seq(0.5,240,by=0.5))) +
geom_line(aes(y = value, col = measure)) +
geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
theme_classic() +
facet_wrap(~measure, scales = "free_y", nrow = 2) +
scale_color_manual(values = c("orange","cornflowerblue")) +
labs(x = "observation")
p_res2 <- residuals(model2, type = "pearson")
d_res2 <- residuals(model2, type = "deviance")
data.frame("pearson" = p_res2, "deviance" = d_res2) %>%
pivot_longer(everything(), names_to = "residual") %>%
ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
geom_point() +
facet_wrap(~residual, nrow = 2) +
theme_classic() +
scale_color_manual(values = c("darkgreen","purple")) +
labs(x = "observation")
data.frame("pearson" = model2.diag$rp, "deviance" = model2.diag$rd) %>%
pivot_longer(everything(), names_to = "residual") %>%
ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
geom_point() +
facet_wrap(~residual, nrow = 2) +
theme_classic() +
scale_color_manual(values = c("darkgoldenrod2","orchid3")) +
labs(x = "observation", y = "standardized values")
#kable(data.frame("pearson" = p_res2, "deviance" = d_res2,"leverage" = model2.diag$h, "cooks" = model2.diag$cook,"p_standard" = model2.diag$rp, "d_standard" = model2.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
# model 3: baseline + selected top genes + immune response genes + age:immunse response genes
model3 <- glm(y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4
+ age_c:CD86 + age_c:DDX58 + age_c:KLRD1 + age_c:TLR4 , data = dd1, family = 'binomial')
summary(model3)
##
## Call:
## glm(formula = y ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c +
## CD86 + DDX58 + KLRD1 + TLR4 + age_c:CD86 + age_c:DDX58 +
## age_c:KLRD1 + age_c:TLR4, family = "binomial", data = dd1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1112 -0.9507 0.3760 0.9613 2.0786
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.291027 0.202136 1.440 0.14993
## MC1R -0.375414 0.295563 -1.270 0.20403
## AKT1 -0.459284 0.228609 -2.009 0.04453 *
## PAK2 0.290848 0.176836 1.645 0.10002
## SOS1 -0.126717 0.159818 -0.793 0.42785
## B2M -0.182347 0.220100 -0.828 0.40740
## sex -0.317183 0.311342 -1.019 0.30832
## age_c -0.029922 0.010791 -2.773 0.00556 **
## CD86 0.295873 0.290758 1.018 0.30887
## DDX58 0.282801 0.210799 1.342 0.17974
## KLRD1 0.053540 0.253813 0.211 0.83293
## TLR4 0.303441 0.233045 1.302 0.19289
## age_c:CD86 -0.006308 0.020972 -0.301 0.76358
## age_c:DDX58 -0.047733 0.019392 -2.461 0.01384 *
## age_c:KLRD1 -0.006409 0.016353 -0.392 0.69510
## age_c:TLR4 0.022713 0.018525 1.226 0.22017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 331.89 on 239 degrees of freedom
## Residual deviance: 266.87 on 224 degrees of freedom
## AIC: 298.87
##
## Number of Fisher Scoring iterations: 5
model3.diag <- glm.diag(model3)
data.frame("cooks" = model3.diag$cook, "leverage" = model3.diag$h) %>%
pivot_longer(everything(), names_to = "measure") %>%
mutate(cutoff = ifelse(measure=="cooks",1,26/240)) %>%
ggplot(aes(x = seq(0.5,240,by=0.5))) +
geom_line(aes(y = value, col = measure)) +
geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
theme_classic() +
facet_wrap(~measure, scales = "free_y", nrow = 2) +
scale_color_manual(values = c("orange","cornflowerblue")) +
labs(x = "observation")
p_res3 <- residuals(model3, type = "pearson")
d_res3 <- residuals(model3, type = "deviance")
data.frame("pearson" = p_res3, "deviance" = d_res3) %>%
pivot_longer(everything(), names_to = "residual") %>%
ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
geom_point() +
facet_wrap(~residual, nrow = 2) +
theme_classic() +
scale_color_manual(values = c("darkgreen","purple")) +
labs(x = "observation")
data.frame("pearson" = model3.diag$rp, "deviance" = model3.diag$rd) %>%
pivot_longer(everything(), names_to = "residual") %>%
ggplot(aes(x = seq(0.5,240,by=0.5), y = value, col = residual)) +
geom_point() +
facet_wrap(~residual, nrow = 2) +
theme_classic() +
scale_color_manual(values = c("darkgoldenrod2","orchid3")) +
labs(x = "observation", y = "standardized values")
#kable(data.frame("pearson" = p_res3, "deviance" = d_res3,"leverage" = model3.diag$h, "cooks" = model3.diag$cook,"p_standard" = model3.diag$rp, "d_standard" = model3.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
library(VGAM)
## Warning: package 'VGAM' was built under R version 4.1.3
## Loading required package: stats4
## Loading required package: splines
##
## Attaching package: 'VGAM'
## The following objects are masked from 'package:boot':
##
## logit, simplex
## The following object is masked from 'package:car':
##
## logit
## The following object is masked from 'package:tidyr':
##
## fill
# model 5: proportional odds model
dd1$CD86 = data$CD86
dd1$DDX58 = data$DDX58
dd1$KLRD1 = data$KLRD1
dd1$TLR4 = data$TLR4
dd1$pathStage = as.factor(data$stage)
model5 <- vglm(pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4, data = dd1, family = cumulative(parallel=T))
## Warning in eval(slot(family, "initialize")): response should be ordinal---see
## ordered()
summary(model5)
##
## Call:
## vglm(formula = pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M +
## sex + age_c + CD86 + DDX58 + KLRD1 + TLR4, family = cumulative(parallel = T),
## data = dd1)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.204631 0.155436 -7.750 9.19e-15 ***
## (Intercept):2 0.288234 0.138098 2.087 0.0369 *
## (Intercept):3 5.690631 1.008460 5.643 1.67e-08 ***
## MC1R -0.318350 0.139875 -2.276 0.0228 *
## AKT1 0.283519 0.124273 2.281 0.0225 *
## PAK2 0.215181 0.104077 2.068 0.0387 *
## SOS1 0.059470 0.098668 0.603 0.5467
## B2M 0.018327 0.127866 0.143 0.8860
## sex -0.179287 0.184433 -0.972 0.3310
## age_c -0.009620 0.005756 -1.671 0.0946 .
## CD86 -0.291645 0.162506 -1.795 0.0727 .
## DDX58 0.245220 0.099962 2.453 0.0142 *
## KLRD1 0.050054 0.136286 0.367 0.7134
## TLR4 -0.030667 0.115717 -0.265 0.7910
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]),
## logitlink(P[Y<=3])
##
## Residual deviance: 504.6351 on 706 degrees of freedom
##
## Log-likelihood: -252.3176 on 706 degrees of freedom
##
## Number of Fisher scoring iterations: 14
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):3'
##
##
## Exponentiated coefficients:
## MC1R AKT1 PAK2 SOS1 B2M sex age_c CD86
## 0.7273479 1.3277943 1.2400857 1.0612743 1.0184961 0.8358656 0.9904259 0.7470341
## DDX58 KLRD1 TLR4
## 1.2779020 1.0513284 0.9697985
model6 <- vglm(pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4
+ age_c:CD86 + age_c:DDX58 + age_c:KLRD1 + age_c:TLR4 , data = dd1, family = cumulative(parallel = T))
## Warning in eval(slot(family, "initialize")): response should be ordinal---see
## ordered()
summary(model6)
##
## Call:
## vglm(formula = pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M +
## sex + age_c + CD86 + DDX58 + KLRD1 + TLR4 + age_c:CD86 +
## age_c:DDX58 + age_c:KLRD1 + age_c:TLR4, family = cumulative(parallel = T),
## data = dd1)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.239885 0.163185 -7.598 3.01e-14 ***
## (Intercept):2 0.325951 0.144582 2.254 0.024168 *
## (Intercept):3 5.822924 1.017363 5.724 1.04e-08 ***
## MC1R -0.386994 0.142384 -2.718 0.006569 **
## AKT1 0.278968 0.128804 2.166 0.030324 *
## PAK2 0.212728 0.107358 1.981 0.047539 *
## SOS1 0.069787 0.102485 0.681 0.495906
## B2M 0.049921 0.137289 0.364 0.716143
## sex -0.221438 0.193277 -1.146 0.251918
## age_c -0.016624 0.006295 -2.641 0.008274 **
## CD86 -0.198693 0.168893 -1.176 0.239420
## DDX58 0.338315 0.115911 2.919 0.003514 **
## KLRD1 -0.058325 0.156900 -0.372 0.710092
## TLR4 -0.090631 0.125003 -0.725 0.468431
## age_c:CD86 0.025759 0.013054 1.973 0.048472 *
## age_c:DDX58 -0.034658 0.010179 -3.405 0.000662 ***
## age_c:KLRD1 -0.027954 0.010426 -2.681 0.007335 **
## age_c:TLR4 -0.015682 0.010922 -1.436 0.151045
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]),
## logitlink(P[Y<=3])
##
## Residual deviance: 489.3675 on 702 degrees of freedom
##
## Log-likelihood: -244.6837 on 702 degrees of freedom
##
## Number of Fisher scoring iterations: 15
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):3'
##
##
## Exponentiated coefficients:
## MC1R AKT1 PAK2 SOS1 B2M sex
## 0.6790952 1.3217656 1.2370476 1.0722794 1.0511879 0.8013654
## age_c CD86 DDX58 KLRD1 TLR4 age_c:CD86
## 0.9835135 0.8198018 1.4025826 0.9433435 0.9133542 1.0260937
## age_c:DDX58 age_c:KLRD1 age_c:TLR4
## 0.9659355 0.9724329 0.9844399